This script analyzes and plots data for Symbiontic Integration 2021 respirometry data for thermal performance curves.

Setup

Set up workspace, set options, and load required packages.

knitr::opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE)
## install packages if you dont already have them in your library
if ("tidyverse" %in% rownames(installed.packages()) == 'FALSE') install.packages('tidyverse') 
if ("car" %in% rownames(installed.packages()) == 'FALSE') install.packages('car') 
if ("scales" %in% rownames(installed.packages()) == 'FALSE') install.packages('scales') 
if ("ggplot2" %in% rownames(installed.packages()) == 'FALSE') install.packages('ggplot2') 
if ("nls.multstart" %in% rownames(installed.packages()) == 'FALSE') install.packages('nls.multstart') 
if ("broom" %in% rownames(installed.packages()) == 'FALSE') install.packages('broom') 
if ("rTPC" %in% rownames(installed.packages()) == 'FALSE') remotes::install_github("padpadpadpad/rTPC")

#load packages
library("ggplot2")
library("tidyverse")
library('car')
library('scales')
library('nls.multstart')
library('broom')
library('rTPC')
library('ggstatsplot')

Data visualization

Load data from LoLinR.

PRdata<-read.csv("Pacu2021/Output/Respiration/TPC/oxygen_Resp_rates.csv") #load data

Format data columns.

#remove all rows of wells that did not have samples or blanks
PRdata<-PRdata[!is.na(PRdata$Type),]

#format columns
PRdata$Temperature<-as.factor(PRdata$Temperature)
PRdata$Temperature<-as.numeric(as.character(PRdata$Temperature))
PRdata$Treatment<-as.factor(PRdata$Treatment)
PRdata$Plate<-as.factor(PRdata$Plate)
PRdata$Run<-as.factor(PRdata$Run)


PRdata <- PRdata %>%
    group_by(Temperature,Treatment)%>%
    mutate(abs=abs(R.nmol.org.min))%>%
  summarise(rate=abs) %>%
  mutate(temp=Temperature)

PRdata$group <- paste0(PRdata$Temperature,"_", PRdata$Treatment)

View and remove outliers.

#PRdata <- PRdata %>% filter(R.nmol.org.min > -0.15) %>% filter(R.nmol.org.min < 0)

#identify outliers by Tremperature and Treatment groups
outlier.plot <- ggbetweenstats(PRdata,group, rate, outlier.tagging = TRUE)
## Warning: Number of labels is greater than default palette color count.
##  Try using another color `palette` (and/or `package`).
outlier.plot
<<<<<<< HEAD

=======

>>>>>>> 2935501d0ff9c1aaecd004f8a449a38a9358f695
ggsave("Pacu2021/Figures/Respiration/TPC/outliersbygroup.pdf", outlier.plot, dpi=300, w=16, h=8, units="in")

q <- c(0.25, 0.75)

Quants <- PRdata %>%
  group_by(Treatment, Temperature) %>%
  summarize(quant25 = quantile(rate, probs = q[1]),
            quant75 = quantile(rate, probs = q[2]),
            IQRbyGroup=IQR(rate))

Quants$group <-paste0(Quants$Temperature,"_", Quants$Treatment)

Quants$upper <-  Quants$quant75+1.5*Quants$IQRbyGroup # Upper Range  
Quants$lower <- Quants$quant25-1.5*Quants$IQRbyGroup # Lower Range

#join outlier cutoffs with data
PRdata <- left_join(PRdata, Quants)

#remove outliers
PRdata <- PRdata %>%
  filter(rate < upper) %>%
  filter(rate > lower)

outlier.removed.plot <- ggbetweenstats(PRdata, group, rate, outlier.tagging = TRUE)
## Warning: Number of labels is greater than default palette color count.
##  Try using another color `palette` (and/or `package`).
outlier.removed.plot
<<<<<<< HEAD

=======

>>>>>>> 2935501d0ff9c1aaecd004f8a449a38a9358f695
ggsave("Pacu2021/Figures/Respiration/TPC/Removed_outliersbygroup.pdf", outlier.removed.plot, dpi=300, w=16, h=8, units="in")

Plotting

Plot data as a dot plot

r_plot<-PRdata %>%
    ggplot(., aes(x = Temperature, y = rate, colour=Treatment)) +
    geom_point(aes(fill=Treatment, group=interaction(Temperature, Treatment)), pch = 21, size=2, alpha=0.5, position = position_jitterdodge(0.3)) + 
    xlab("Temperature") + 
    scale_fill_manual(name="Treatment", values=c("blue","red"))+
    scale_color_manual(name="Treatment", values=c("blue","red"))+
    ylab(expression(bold(paste("R (nmol ", O[2], " larva"^-1, "min"^-1, ")")))) +
    theme_classic() + 
    theme(
      legend.position="none",
      axis.title=element_text(face="bold", size=16),
      axis.text=element_text(size=12, color="black"), 
      legend.title=element_text(face="bold", size=14), 
      legend.text=element_text(size=12)
      ); r_plot
<<<<<<< HEAD

=======

>>>>>>> 2935501d0ff9c1aaecd004f8a449a38a9358f695
ggsave("Pacu2021/Figures/Respiration/TPC/Respiration_Scatter.pdf", r_plot, dpi=300, w=8, h=8, units="in")

Filtering

PRdata <- PRdata %>%
    filter(!Temperature==8) %>%
    filter(!Temperature==18)%>%
    filter(!Temperature==34)%>%
    filter(!Temperature==36)%>%
    filter(!Temperature==38)%>%
    filter(!Temperature==40)

r_plot<-PRdata %>%
    ggplot(., aes(x = Temperature, y = rate, colour=Treatment)) +
    geom_point(aes(fill=Treatment, group=interaction(Temperature, Treatment)), pch = 21, size=2, alpha=0.5, position = position_jitterdodge(0.3)) + 
   #geom_smooth(method = "lm", se=FALSE, formula=y ~ poly(x, 3, raw=TRUE), size=3)+
    xlab("Temperature") + 
    scale_fill_manual(name="Treatment", values=c("blue","red"))+
    scale_color_manual(name="Treatment", values=c("blue","red"))+
    ylab(expression(bold(paste("R (nmol ", O[2], " larva"^-1, "min"^-1, ")")))) +
    theme_classic() + 
    theme(
      legend.position="none",
      axis.title=element_text(face="bold", size=16),
      axis.text=element_text(size=12, color="black"), 
      legend.title=element_text(face="bold", size=14), 
      legend.text=element_text(size=12)
      ); r_plot
<<<<<<< HEAD

=======

>>>>>>> 2935501d0ff9c1aaecd004f8a449a38a9358f695
ggsave("Pacu2021/Figures/Respiration/TPC/Respiration_Scatter_Filtered.pdf", r_plot, dpi=300, w=8, h=8, units="in")

View and remove outliers.

#identify outliers by Tremperature and Treatment groups
outlier.plot <- ggbetweenstats(PRdata,group, rate, outlier.tagging = TRUE)
## Warning: Number of labels is greater than default palette color count.
##  Try using another color `palette` (and/or `package`).
outlier.plot
<<<<<<< HEAD

=======

>>>>>>> 2935501d0ff9c1aaecd004f8a449a38a9358f695
Quants <- PRdata %>%
  group_by(Treatment, Temperature) %>%
  summarize(quant25 = quantile(rate, probs = q[1]),
            quant75 = quantile(rate, probs = q[2]),
            IQRbyGroup=IQR(rate))

Quants$group <-paste0(Quants$Temperature,"_", Quants$Treatment)

Quants$upper <-  Quants$quant75+1.5*Quants$IQRbyGroup # Upper Range
Quants$lower <- Quants$quant25-1.5*Quants$IQRbyGroup # Lower Range

#join outlier cutoffs with data
PRdata <- left_join(PRdata, Quants)

#remove outliers
PRdata <- PRdata %>%
  filter(rate < upper) %>%
  filter(rate > lower)

<<<<<<< HEAD
PRdata <- PRdata %>% filter(rate > 0.0015)
=======
>>>>>>> 2935501d0ff9c1aaecd004f8a449a38a9358f695
PRdata <- PRdata %>% filter(rate < 0.14)
PRdata <- PRdata[-c(96,106),]

outlier.removed.plot <- ggbetweenstats(PRdata, group, rate, outlier.tagging = TRUE)
## Warning: Number of labels is greater than default palette color count.
##  Try using another color `palette` (and/or `package`).
outlier.removed.plot
<<<<<<< HEAD

=======

>>>>>>> 2935501d0ff9c1aaecd004f8a449a38a9358f695
ggsave("Pacu2021/Figures/Respiration/TPC/Removed_outliersbygroup2.pdf", outlier.removed.plot, dpi=300, w=16, h=8, units="in")

Analysis

TPC ftting

#Padifeld et al rTPC and nls.multstart: A new pipeline to fit thermal performance curves in r ### https://besjournals.onlinelibrary.wiley.com/doi/full/10.1111/2041-210X.13585

# choose model
get_model_names()
##  [1] "beta_2012"             "boatman_2017"          "briere2_1999"         
##  [4] "delong_2017"           "flinn_1991"            "gaussian_1987"        
##  [7] "hinshelwood_1947"      "joehnk_2008"           "johnsonlewin_1946"    
## [10] "kamykowski_1985"       "lactin2_1995"          "lrf_1991"             
## [13] "modifiedgaussian_2006" "oneill_1972"           "pawar_2018"           
## [16] "quadratic_2008"        "ratkowsky_1983"        "rezende_2019"         
## [19] "sharpeschoolfull_1981" "sharpeschoolhigh_1981" "sharpeschoollow_1981" 
## [22] "spain_1982"            "thomas_2012"           "thomas_2017"          
## [25] "weibull_1995"
#sharpeschoolhigh_1981

# get start vals
start_vals <- get_start_vals(PRdata$temp,PRdata$rate, model_name = 'sharpeschoolhigh_1981')

# get limits
low_lims <- get_lower_lims(PRdata$temp,PRdata$rate, model_name = 'sharpeschoolhigh_1981')
upper_lims <- get_upper_lims(PRdata$temp,PRdata$rate, model_name = 'sharpeschoolhigh_1981')

start_vals
<<<<<<< HEAD
##     r_tref          e         eh         th 
##  0.0456096  0.7027915  2.1208907 30.9411765
=======
##      r_tref           e          eh          th 
##  0.04434088  0.76808284  2.42820170 30.94117647
>>>>>>> 2935501d0ff9c1aaecd004f8a449a38a9358f695
low_lims
## r_tref      e     eh     th 
##      0      0      0      1
upper_lims
##     r_tref          e         eh         th 
##  0.1232836 10.0000000 20.0000000 32.0000000
# AMBIENT CURVE FIT
d.amb <- PRdata %>% 
  filter(Treatment=="Ambient")

# # get start vals
# amb.start_vals <- get_start_vals(d.amb$temp,d.amb$rate, model_name = 'sharpeschoolhigh_1981')
# amb.start_vals
# 
# # get limits
# amb.low_lims <- get_lower_lims(d.amb$temp,d.amb$rate, model_name = 'sharpeschoolhigh_1981')
# amb.low_lims
# amb.upper_lims <- get_upper_lims(d.amb$temp,d.amb$rate, model_name = 'sharpeschoolhigh_1981')
# amb.start_vals
# amb.low_lims
# amb.upper_lims

amb.fit <- nls_multstart(rate~sharpeschoolhigh_1981(temp = temp, r_tref,e,eh,th, tref = 28),
                                                     data = d.amb,
                                                     iter = 500,
                                                     start_lower = start_vals - 1,
                                                     start_upper = start_vals + 1,
                                                     lower = low_lims,
                                                     upper = upper_lims,
                                                     supp_errors = 'Y')

amb.fit
## Nonlinear regression model
##   model: rate ~ sharpeschoolhigh_1981(temp = temp, r_tref, e, eh, th,     tref = 28)
##    data: data
##   r_tref        e       eh       th 
<<<<<<< HEAD
##  0.06369  0.77099  8.81372 31.66205 
##  residual sum-of-squares: 0.03744
## 
## Number of iterations to convergence: 25 
=======
##  0.06287  0.81631 10.05527 31.49888 
##  residual sum-of-squares: 0.03564
## 
## Number of iterations to convergence: 21 
>>>>>>> 2935501d0ff9c1aaecd004f8a449a38a9358f695
## Achieved convergence tolerance: 1.49e-08
# predict new data
amb_new_data <- data.frame(temp = seq(min(d.amb$temp), max(d.amb$temp), 0.5))
amb.preds <- augment(amb.fit, newdata = amb_new_data)

amb.TCP.res <- calc_params(amb.fit) %>%
  mutate_all(round, 2)   # round 

amb.TCP.res 
<<<<<<< HEAD
##   rmax  topt ctmin ctmax   e   eh q10 thermal_safety_margin thermal_tolerance
## 1 0.07 29.55  0.81 34.96 0.7 2.31 2.5                  5.41             34.15
##   breadth skewness
## 1    4.99    -1.61
=======
##   rmax  topt ctmin ctmax    e  eh  q10 thermal_safety_margin thermal_tolerance
## 1 0.07 29.58  2.34 34.35 0.77 2.9 2.71                  4.77             32.01
##   breadth skewness
## 1    4.58    -2.14
>>>>>>> 2935501d0ff9c1aaecd004f8a449a38a9358f695

BOOTSTRAPPING AMBIENT CURVE FIT

Resampling the original data with replacement

# refit model using nlsLM
amb.fit_nlsLM <- minpack.lm::nlsLM(rate~sharpeschoolhigh_1981(temp = temp, r_tref,e,eh,th, tref = 28),
                        data = d.amb,
                        start = coef(amb.fit),
                        lower = low_lims,
                        upper = upper_lims,
                        weights = rep(1, times = nrow(d.amb)))

# bootstrap using case resampling
amb.boot1 <- Boot(amb.fit_nlsLM, method = 'case')

# look at the data
head(amb.boot1$t)
##          r_tref         e        eh       th
<<<<<<< HEAD
## [1,] 0.05622679 1.0115372  3.274674 32.00000
## [2,] 0.06621550 0.6028136 12.494659 31.84514
## [3,] 0.11625218 1.5155506  4.932409 29.19690
## [4,] 0.12328360 1.6623770  3.525770 28.08559
## [5,] 0.12328360 1.5686641  2.828922 27.31338
## [6,] 0.10076807 1.2078562  3.945774 29.13973
======= ## [1,] 0.12328360 1.8657450 4.027641 28.04144 ## [2,] 0.09298738 2.4680678 7.536008 29.25540 ## [3,] 0.06189850 0.8984202 20.000000 31.68621 ## [4,] 0.12328360 2.0854617 5.215647 28.79464 ## [5,] 0.07015725 1.0264364 20.000000 31.58002 ## [6,] 0.07430899 1.1065485 9.812652 31.17470 >>>>>>> 2935501d0ff9c1aaecd004f8a449a38a9358f695
# create predictions of each bootstrapped model
amb.boot1_preds <- amb.boot1$t %>%
  as.data.frame() %>%
  drop_na() %>%
  mutate(iter = 1:n()) %>%
  group_by_all() %>%
  do(data.frame(temp = seq(min(d.amb$temp), max(d.amb$temp), length.out = 100))) %>%
  ungroup() %>%
  mutate(pred = sharpeschoolhigh_1981(temp, r_tref, e, eh, th, tref = 28))

# calculate bootstrapped confidence intervals
amb.boot1_conf_preds <- group_by(amb.boot1_preds, temp) %>%
  summarise(conf_lower = quantile(pred, 0.025),
            conf_upper = quantile(pred, 0.975)) %>%
  ungroup()

# plot bootstrapped CIs
amb.CI.plot <- ggplot() +
  geom_line(aes(temp, .fitted), amb.preds, col = 'blue') +
  geom_ribbon(aes(temp, ymin = conf_lower, ymax = conf_upper), amb.boot1_conf_preds, fill = 'blue', alpha = 0.3) +
  geom_point(aes(temp, rate), d.amb, size = 2, alpha = 0.5,col = 'blue') +
  theme_bw(base_size = 12) +
  labs(x = 'Temperature (ºC)',
       y = 'Respiration Rate (nmol O2/larva/min')
amb.CI.plot 
<<<<<<< HEAD

=======

>>>>>>> 2935501d0ff9c1aaecd004f8a449a38a9358f695
# HIGH CURVE FIT
d.high <- PRdata %>% 
  filter(Treatment=="High")

# choose model
mod <- 'sharpschoolhigh_1981'

# # get start vals
# high.start_vals <- get_start_vals(d.high$temp,d.high$rate, model_name = 'sharpeschoolhigh_1981')
# high.start_vals
# 
# # get limits
# high.low_lims <- get_lower_lims(d.high$temp,d.high$rate, model_name = 'sharpeschoolhigh_1981')
# high.low_lims
# high.upper_lims <- get_upper_lims(d.high$temp,d.high$rate, model_name = 'sharpeschoolhigh_1981')
# high.start_vals
# high.low_lims
# high.upper_lims

high.fit <- nls_multstart(rate~sharpeschoolhigh_1981(temp = temp, r_tref,e,eh,th, tref = 28),
                                                     data = d.high,
                                                     iter = 500,
                                                     start_lower = start_vals - 1,
                                                     start_upper = start_vals + 1,
                                                     lower = low_lims,
                                                     upper = upper_lims,
                                                     supp_errors = 'Y')

high.fit
## Nonlinear regression model
##   model: rate ~ sharpeschoolhigh_1981(temp = temp, r_tref, e, eh, th,     tref = 28)
##    data: data
##  r_tref       e      eh      th 
##  0.1233  1.4038  2.4282 26.4288 
##  residual sum-of-squares: 0.03235
## 
## Number of iterations to convergence: 9 
## Achieved convergence tolerance: 1.49e-08
high_new_data <- data.frame(temp = seq(min(d.high$temp), max(d.high$temp), 0.5))
high.preds <- augment(high.fit, newdata = high_new_data)


high.TCP.res <- calc_params(high.fit) %>%
  mutate_all(round, 2)   # round for easy viewing

high.TCP.res 
##   rmax  topt ctmin ctmax    e   eh  q10 thermal_safety_margin thermal_tolerance
## 1 0.05 27.44  7.45 57.05 0.53 0.47 2.01                 29.62              49.6
##   breadth skewness
## 1    8.87     0.06

BOOTSTRAPPING HIGH CURVE FIT

Resampling the original data with replacement

# refit model using nlsLM
high.fit_nlsLM <- minpack.lm::nlsLM(rate~sharpeschoolhigh_1981(temp = temp, r_tref,e,eh,th, tref = 28),
                        data = d.high,
                        start = coef(high.fit),
                        lower = low_lims,
                        upper = upper_lims,
                        weights = rep(1, times = nrow(d.high)))

# bootstrap using case resampling
<<<<<<< HEAD
high.boot1 <- Boot(high.fit_nlsLM, method = 'case')

# look at the data
head(high.boot1$t)
##         r_tref         e        eh       th
## [1,] 0.1232836 1.0875773 1.3645045 24.82696
## [2,] 0.1232836 1.0663405 1.7193422 25.18904
## [3,] 0.1232836 1.5667784 3.4832088 27.14641
## [4,] 0.1232836 1.2804056 2.9616468 26.35147
## [5,] 0.1232836 1.5242106 3.1892471 27.27148
## [6,] 0.1232836 0.0798746 0.1515831  1.00000
======= high.boot1 <- Boot(high.fit_nlsLM, method = 'case')
## 
##  Number of bootstraps was 997 out of 999 attempted
# look at the data
head(high.boot1$t)
##          r_tref         e       eh       th
## [1,] 0.09466776 1.1017541 3.245382 28.74252
## [2,] 0.12328360 1.5033158 2.816571 26.55534
## [3,] 0.12328360 1.3592890 2.435695 26.91997
## [4,] 0.12328360 1.6925379 2.945170 26.25989
## [5,] 0.07216997 0.8133763 4.153968 29.78028
## [6,] 0.04920720 0.3859032 4.080503 32.00000
>>>>>>> 2935501d0ff9c1aaecd004f8a449a38a9358f695
# create predictions of each bootstrapped model
high.boot1_preds <- high.boot1$t %>%
  as.data.frame() %>%
  drop_na() %>%
  mutate(iter = 1:n()) %>%
  group_by_all() %>%
  do(data.frame(temp = seq(min(d.high$temp), max(d.high$temp), length.out = 100))) %>%
  ungroup() %>%
  mutate(pred = sharpeschoolhigh_1981(temp, r_tref, e, eh, th, tref = 28))

# calculate bootstrapped confidence intervals
high.boot1_conf_preds <- group_by(high.boot1_preds, temp) %>%
  summarise(conf_lower = quantile(pred, 0.025),
            conf_upper = quantile(pred, 0.975)) %>%
  ungroup()

# plot bootstrapped CIs
high.CI.plot <- ggplot() +
  geom_line(aes(temp, .fitted), high.preds, col = 'red') +
  geom_ribbon(aes(temp, ymin = conf_lower, ymax = conf_upper), high.boot1_conf_preds, fill = 'red', alpha = 0.3) +
  geom_point(aes(temp, rate), d.high, size = 2, alpha = 0.5,col = 'red') +
  theme_bw(base_size = 12) +
  labs(x = 'Temperature (ºC)',
       y = 'Respiration Rate (nmol O2/larva/min')
high.CI.plot
<<<<<<< HEAD

=======

>>>>>>> 2935501d0ff9c1aaecd004f8a449a38a9358f695
# plot data and model fit
TPC.plot <- ggplot() +
   geom_point(aes(temp, rate), d.amb, size = 2, alpha = 0.5,col = 'blue') +
   geom_point(aes(temp, rate), d.high, size = 2, alpha = 0.5,col = 'red') +
   geom_line(aes(temp, .fitted), amb.preds, col = 'blue', size=2) +
   geom_line(aes(temp, .fitted), high.preds, col = 'red', size=2) +
   geom_ribbon(aes(temp, ymin = conf_lower, ymax = conf_upper), high.boot1_conf_preds, fill = 'red', alpha = 0.3) +
   geom_ribbon(aes(temp, ymin = conf_lower, ymax = conf_upper), amb.boot1_conf_preds, fill = 'blue', alpha = 0.3) +
   theme_bw(base_size = 12) +
   labs(x = 'Temperature (ºC)',
        y = 'Metabolic rate')

TPC.plot 
<<<<<<< HEAD

ggsave("Pacu2021/Figures/Respiration/TPC/TPC.pdf", TPC.plot, dpi=300, w=8, h=8, units="in")
=======

ggsave("Pacu2021/Figures/Respiration/TPC/TPC_SharpSchool.pdf", TPC.plot, dpi=300, w=8, h=8, units="in")
>>>>>>> 2935501d0ff9c1aaecd004f8a449a38a9358f695

#confidence intervals of TPC parameters

broom::tidy(amb.fit_nlsLM)
## # A tibble: 4 x 5
##   term   estimate std.error statistic  p.value
##   <chr>     <dbl>     <dbl>     <dbl>    <dbl>
<<<<<<< HEAD
## 1 r_tref   0.0637    0.0110      5.78 5.07e- 7
## 2 e        0.771     0.412       1.87 6.76e- 2
## 3 eh       8.81      7.46        1.18 2.43e- 1
## 4 th      31.7       0.850      37.2  1.34e-37
======= ## 1 r_tref 0.0629 0.00968 6.50 3.67e- 8 ## 2 e 0.816 0.385 2.12 3.89e- 2 ## 3 eh 10.1 7.32 1.37 1.75e- 1 ## 4 th 31.5 0.755 41.7 1.56e-40 >>>>>>> 2935501d0ff9c1aaecd004f8a449a38a9358f695
broom::tidy(high.fit_nlsLM)
## # A tibble: 4 x 5
##   term   estimate std.error statistic p.value
##   <chr>     <dbl>     <dbl>     <dbl>   <dbl>
## 1 r_tref    0.123     0.497     0.248  0.805 
## 2 e         1.40      3.76      0.373  0.710 
## 3 eh        2.43      1.14      2.12   0.0384
## 4 th       26.4      20.7       1.28   0.208
#AMBIENT
amb.extra_params <- calc_params(amb.fit_nlsLM) %>%
  pivot_longer(everything(), names_to =  'param', values_to = 'estimate')

amb.ci_extra_params <- Boot(amb.fit_nlsLM, f = function(x){unlist(calc_params(x))}, labels = names(calc_params(amb.fit_nlsLM)), R = 200, method = 'case') %>%
  confint(., method = 'bca') %>%
  as.data.frame() %>%
  rename(conf_lower = 1, conf_upper = 2) %>%
  rownames_to_column(., var = 'param') %>%
  mutate(method = 'case bootstrap')
## 
<<<<<<< HEAD
##  Number of bootstraps was 83 out of 200 attempted
======= ## Number of bootstraps was 68 out of 200 attempted >>>>>>> 2935501d0ff9c1aaecd004f8a449a38a9358f695
amb.ci_extra_params <- left_join(amb.ci_extra_params, amb.extra_params)
amb.ci_extra_params$Treatment <- "Ambient"
# Joining, by = "param"


#HIGH
high.extra_params <- calc_params(high.fit_nlsLM) %>%
  pivot_longer(everything(), names_to =  'param', values_to = 'estimate')

high.ci_extra_params <- Boot(high.fit_nlsLM, f = function(x){unlist(calc_params(x))}, labels = names(calc_params(high.fit_nlsLM)), R = 200, method = 'case') %>%
  confint(., method = 'bca') %>%
  as.data.frame() %>%
  rename(conf_lower = 1, conf_upper = 2) %>%
  rownames_to_column(., var = 'param') %>%
  mutate(method = 'case bootstrap')
## 
<<<<<<< HEAD
##  Number of bootstraps was 139 out of 200 attempted
======= ## Number of bootstraps was 151 out of 200 attempted
high.ci_extra_params <- left_join(high.ci_extra_params, high.extra_params)
high.ci_extra_params$Treatment <- "High"

#Join Ambient and High estimates and CIs

All_params <- rbind(amb.ci_extra_params, high.ci_extra_params)
All_params <- All_params %>% 
 mutate_if(is.numeric, round, 2)

#Plot
estimate.plots <- ggplot(All_params, aes(Treatment, estimate, color=Treatment)) +
  geom_point(size = 2) +
  scale_color_manual(name="Treatment", values=c("blue","red"))+
  geom_linerange(aes(ymin = conf_lower, ymax = conf_upper)) +
  theme_bw() +
  facet_wrap(~param, scales = 'free_y') +
  scale_x_discrete('')

estimate.plots

ggsave("Pacu2021/Figures/Respiration/TPC/TPC_estimates_SharpSchool.pdf", estimate.plots, dpi=300, w=8, h=8, units="in")

#Gaussian

#AMBIENT
amb.gaus.fit <- nls_multstart(rate~gaussian_1987(temp = temp, rmax, topt, a),
                        data = d.amb,
                        iter = c(4,4,4),
                        start_lower = get_start_vals(d.amb$temp, d.amb$rate, model_name = 'gaussian_1987') - 1,
                        start_upper = get_start_vals(d.amb$temp, d.amb$rate, model_name = 'gaussian_1987') + 1,
                        lower = get_lower_lims(d.amb$temp, d.amb$rate, model_name = 'gaussian_1987'),
                        upper = get_upper_lims(d.amb$temp, d.amb$rate, model_name = 'gaussian_1987'),
                        supp_errors = 'Y',
                        convergence_count = FALSE)


amb.gaus.fit 
## Nonlinear regression model
##   model: rate ~ gaussian_1987(temp = temp, rmax, topt, a)
##    data: data
##     rmax     topt        a 
##  0.06347 27.83799  4.38800 
##  residual sum-of-squares: 0.03694
## 
## Number of iterations to convergence: 11 
## Achieved convergence tolerance: 1.49e-08
# predict new data
amb_new_data <- data.frame(temp = seq(min(d.amb$temp), max(d.amb$temp), 0.5))
amb.preds <- augment(amb.gaus.fit, newdata = amb_new_data)

amb.TCP.res <- calc_params(amb.gaus.fit) %>%
  mutate_all(round, 2)   # round for easy viewing
amb.TCP.res 
##   rmax  topt ctmin ctmax    e   eh q10 thermal_safety_margin thermal_tolerance
## 1 0.06 27.84  17.1 38.58 1.37 0.73  NA                 10.74             21.48
##   breadth skewness
## 1    5.87     0.65
amb.gaus.fit_nlsLM <- minpack.lm::nlsLM(rate~gaussian_1987(temp = temp, rmax, topt, a),
                        data = d.amb,
                        start = coef(amb.gaus.fit),
                        lower = get_lower_lims(d.amb$temp, d.amb$rate, model_name = 'gaussian_1987'),
                        upper = get_upper_lims(d.amb$temp, d.amb$rate, model_name = 'gaussian_1987'),
                        weights = rep(1, times = nrow(d.amb)))

# bootstrap using case resampling
amb.boot1 <- Boot(amb.gaus.fit_nlsLM, method = 'case')

# look at the data
head(amb.boot1$t)
##            rmax     topt        a
## [1,] 0.06197087 27.13551 4.854462
## [2,] 0.05900555 27.60181 4.673603
## [3,] 0.05973728 28.41771 4.843638
## [4,] 0.06388289 27.78184 4.564628
## [5,] 0.05356135 27.06235 5.075532
## [6,] 0.07453052 28.19593 3.481052
# create predictions of each bootstrapped model
amb.boot1_preds <- amb.boot1$t %>%
  as.data.frame() %>%
  drop_na() %>%
  mutate(iter = 1:n()) %>%
  group_by_all() %>%
  do(data.frame(temp = seq(min(d.amb$temp), max(d.amb$temp), length.out = 100))) %>%
  ungroup() %>%
  mutate(pred = gaussian_1987(temp = temp, rmax, topt, a))

# calculate bootstrapped confidence intervals
amb.boot1_conf_preds <- group_by(amb.boot1_preds, temp) %>%
  summarise(conf_lower = quantile(pred, 0.025),
            conf_upper = quantile(pred, 0.975)) %>%
  ungroup()

# plot bootstrapped CIs
amb.CI.plot <- ggplot() +
  geom_line(aes(temp, .fitted), amb.preds, col = 'blue') +
  geom_ribbon(aes(temp, ymin = conf_lower, ymax = conf_upper), amb.boot1_conf_preds, fill = 'blue', alpha = 0.3) +
  geom_point(aes(temp, rate), d.amb, size = 2, alpha = 0.5,col = 'blue') +
  theme_bw(base_size = 12) +
  labs(x = 'Temperature (ºC)',
       y = 'Respiration Rate (nmol O2/larva/min')
amb.CI.plot 

# HIGH
high.gaus.fit <- nls_multstart(rate~gaussian_1987(temp = temp, rmax, topt, a),
                        data = d.high,
                        iter = c(4,4,4),
                        start_lower = get_start_vals(d.high$temp, d.high$rate, model_name = 'gaussian_1987') - 1,
                        start_upper = get_start_vals(d.high$temp, d.high$rate, model_name = 'gaussian_1987') + 1,
                        lower = get_lower_lims(d.high$temp, d.high$rate, model_name = 'gaussian_1987'),
                        upper = get_upper_lims(d.high$temp, d.high$rate, model_name = 'gaussian_1987'),
                        supp_errors = 'Y',
                        convergence_count = FALSE)


high.gaus.fit 
## Nonlinear regression model
##   model: rate ~ gaussian_1987(temp = temp, rmax, topt, a)
##    data: data
##     rmax     topt        a 
##  0.04746 27.07125  5.77867 
##  residual sum-of-squares: 0.0274
## 
## Number of iterations to convergence: 9 
## Achieved convergence tolerance: 1.49e-08
# predict new data
high_new_data <- data.frame(temp = seq(min(d.high$temp), max(d.high$temp), 0.5))
high.preds <- augment(high.gaus.fit, newdata = high_new_data)

high.TCP.res <- calc_params(high.gaus.fit) %>%
  mutate_all(round, 2)   # round for easy viewing
high.TCP.res 
##   rmax  topt ctmin ctmax    e   eh  q10 thermal_safety_margin thermal_tolerance
## 1 0.05 27.07 12.93 41.22 0.53 0.79 2.01                 14.15             28.29
##   breadth skewness
## 1    7.72    -0.26
high.gaus.fit_nlsLM <- minpack.lm::nlsLM(rate~gaussian_1987(temp = temp, rmax, topt, a),
                        data = d.high,
                        start = coef(high.gaus.fit),
                        lower = get_lower_lims(d.high$temp, d.high$rate, model_name = 'gaussian_1987'),
                        upper = get_upper_lims(d.high$temp, d.high$rate, model_name = 'gaussian_1987'),
                        weights = rep(1, times = nrow(d.high)))

# bootstrap using case resampling
high.boot1 <- Boot(high.gaus.fit_nlsLM, method = 'case')

# look at the data
head(high.boot1$t)
##            rmax     topt        a
## [1,] 0.05927776 26.86467 4.041476
## [2,] 0.05146475 26.79550 4.483799
## [3,] 0.04765019 30.66779 8.565635
## [4,] 0.04808687 28.83889 8.571844
## [5,] 0.04426714 25.74916 8.505839
## [6,] 0.04689024 28.19428 6.521274
# create predictions of each bootstrapped model
high.boot1_preds <- high.boot1$t %>%
  as.data.frame() %>%
  drop_na() %>%
  mutate(iter = 1:n()) %>%
  group_by_all() %>%
  do(data.frame(temp = seq(min(d.high$temp), max(d.high$temp), length.out = 100))) %>%
  ungroup() %>%
  mutate(pred = gaussian_1987(temp = temp, rmax, topt, a))

# calculate bootstrapped confidence intervals
high.boot1_conf_preds <- group_by(high.boot1_preds, temp) %>%
  summarise(conf_lower = quantile(pred, 0.025),
            conf_upper = quantile(pred, 0.975)) %>%
  ungroup()

# plot bootstrapped CIs
high.CI.plot <- ggplot() +
  geom_line(aes(temp, .fitted), high.preds, col = 'red') +
  geom_ribbon(aes(temp, ymin = conf_lower, ymax = conf_upper), high.boot1_conf_preds, fill = 'red', alpha = 0.3) +
  geom_point(aes(temp, rate), d.high, size = 2, alpha = 0.5,col = 'red') +
  theme_bw(base_size = 12) +
  labs(x = 'Temperature (ºC)',
       y = 'Respiration Rate (nmol O2/larva/min')
high.CI.plot 

# plot data and model fit
TPC.plot <- ggplot() +
   geom_point(aes(temp, rate), d.amb, size = 2, alpha = 0.5,col = 'blue') +
   geom_point(aes(temp, rate), d.high, size = 2, alpha = 0.5,col = 'red') +
   geom_line(aes(temp, .fitted), amb.preds, col = 'blue', size=2) +
   geom_line(aes(temp, .fitted), high.preds, col = 'red', size=2) +
   geom_ribbon(aes(temp, ymin = conf_lower, ymax = conf_upper), high.boot1_conf_preds, fill = 'red', alpha = 0.3) +
   geom_ribbon(aes(temp, ymin = conf_lower, ymax = conf_upper), amb.boot1_conf_preds, fill = 'blue', alpha = 0.3) +
   theme_bw(base_size = 12) +
   labs(x = 'Temperature (ºC)',
        y = 'Metabolic rate')

TPC.plot 

ggsave("Pacu2021/Figures/Respiration/TPC/TPC_Gaussian.pdf", TPC.plot, dpi=300, w=8, h=8, units="in")


amb.TCP.res 
##   rmax  topt ctmin ctmax    e   eh q10 thermal_safety_margin thermal_tolerance
## 1 0.06 27.84  17.1 38.58 1.37 0.73  NA                 10.74             21.48
##   breadth skewness
## 1    5.87     0.65
high.TCP.res 
##   rmax  topt ctmin ctmax    e   eh  q10 thermal_safety_margin thermal_tolerance
## 1 0.05 27.07 12.93 41.22 0.53 0.79 2.01                 14.15             28.29
##   breadth skewness
## 1    7.72    -0.26

#confidence intervals of Gaussian TPC parameters

broom::tidy(amb.gaus.fit_nlsLM)
## # A tibble: 3 x 5
##   term  estimate std.error statistic  p.value
##   <chr>    <dbl>     <dbl>     <dbl>    <dbl>
## 1 rmax    0.0635   0.00604     10.5  2.38e-14
## 2 topt   27.8      0.570       48.9  1.67e-44
## 3 a       4.39     0.778        5.64 7.48e- 7
broom::tidy(high.gaus.fit_nlsLM)
## # A tibble: 3 x 5
##   term  estimate std.error statistic  p.value
##   <chr>    <dbl>     <dbl>     <dbl>    <dbl>
## 1 rmax    0.0475   0.00492      9.65 2.91e-13
## 2 topt   27.1      0.847       32.0  2.56e-36
## 3 a       5.78     1.60         3.61 6.71e- 4
#AMBIENT
amb.extra_params <- calc_params(amb.gaus.fit_nlsLM) %>%
  pivot_longer(everything(), names_to =  'param', values_to = 'estimate')

amb.ci_extra_params <- Boot(amb.gaus.fit_nlsLM, f = function(x){unlist(calc_params(x))}, labels = names(calc_params(amb.gaus.fit_nlsLM)), R = 200, method = 'case') %>%
  confint(., method = 'bca') %>%
  as.data.frame() %>%
  rename(conf_lower = 1, conf_upper = 2) %>%
  rownames_to_column(., var = 'param') %>%
  mutate(method = 'case bootstrap')
## 
##  Number of bootstraps was 109 out of 200 attempted
amb.ci_extra_params <- left_join(amb.ci_extra_params, amb.extra_params)
amb.ci_extra_params$Treatment <- "Ambient"
# Joining, by = "param"


#HIGH
high.extra_params <- calc_params(high.gaus.fit_nlsLM) %>%
  pivot_longer(everything(), names_to =  'param', values_to = 'estimate')

high.ci_extra_params <- Boot(high.gaus.fit_nlsLM, f = function(x){unlist(calc_params(x))}, labels = names(calc_params(high.gaus.fit_nlsLM)), R = 200, method = 'case') %>%
  confint(., method = 'bca') %>%
  as.data.frame() %>%
  rename(conf_lower = 1, conf_upper = 2) %>%
  rownames_to_column(., var = 'param') %>%
  mutate(method = 'case bootstrap')
## 
##  Number of bootstraps was 159 out of 200 attempted
>>>>>>> 2935501d0ff9c1aaecd004f8a449a38a9358f695
high.ci_extra_params <- left_join(high.ci_extra_params, high.extra_params)
high.ci_extra_params$Treatment <- "High"

#Join Ambient and High estimates and CIs

All_params <- rbind(amb.ci_extra_params, high.ci_extra_params)
All_params <- All_params %>% 
 mutate_if(is.numeric, round, 2)

#Plot
estimate.plots <- ggplot(All_params, aes(Treatment, estimate, color=Treatment)) +
  geom_point(size = 2) +
  scale_color_manual(name="Treatment", values=c("blue","red"))+
  geom_linerange(aes(ymin = conf_lower, ymax = conf_upper)) +
  theme_bw() +
  facet_wrap(~param, scales = 'free_y') +
  scale_x_discrete('')

estimate.plots
<<<<<<< HEAD

ggsave("Pacu2021/Figures/Respiration/TPC/TPC_estimates.pdf", estimate.plots, dpi=300, w=8, h=8, units="in")
=======

ggsave("Pacu2021/Figures/Respiration/TPC/TPC_estimates_Gaussian.pdf", estimate.plots, dpi=300, w=8, h=8, units="in")
>>>>>>> 2935501d0ff9c1aaecd004f8a449a38a9358f695